##############
##############
###Figure 2###
##############
##############
# Load necessary packages
library(dplyr)
library(ggplot2)
library(purrr)
library(stringr)
library(tidyverse)

###############################################
###Pull and reorganize the data for plotting### 
###############################################

#Extract simulation values
all_sims <- list(sim_results_p, sim_results_p2,sim_results_env, sim_results_env2, sim_results_cy, sim_results_cy2, sim_results_cs, sim_results_cs2, sim_results_hs, sim_results_hs2)
#Create area labels
area_labels <- c("Agricultural", "Forest")
# Summary function
compute_summary <- function(values) {
  mean_val <- mean(values)
  lower_ci <- quantile(values, 0.025)
  upper_ci <- quantile(values, 0.975)
  return(c(Mean = mean_val, Lower = lower_ci, Upper = upper_ci))
}

#Draw the results and calculate 95% CIs
summary_list <- mapply(function(sim_results, area_label) {
  ACME <- unlist(lapply(sim_results, function(x) x$d.avg.sims))
  
  summary_matrix <- t(sapply(list(ACME), compute_summary))
  colnames(summary_matrix) <- c("Mean", "Lower", "Upper")
  
  df <- data.frame(
    Area = area_label,
    Effect = c("ACME"),
    summary_matrix
  )
  return(df)
}, all_sims , area_labels, SIMPLIFY = FALSE)

#Reorganize the resulting data file
summary_all <- do.call(rbind, summary_list)
summary_all$Test <- factor(summary_all$Effect, levels =  c("ACME")) # c("Placebo","CY", "CSE", "HR"))
summary_all$Area <- factor(summary_all$Area, levels = c("Agricultural", "Forest"))
summary_all$row_order <- with(summary_all, interaction(Test, Area, drop = TRUE))
summary_all$row_order <-  c(
  "Placebo.Agricultural", "Placebo.Forest",
  "ENV.Agricultural", "ENV.Forest",
  "CY.Agricultural", "CY.Forest",
  "CSE.Agricultural", "CSE.Forest", 
  "HR.Agricultural", "HR.Forest"
)

#Convert to percents
summary_percent <- summary_all
#Make sure it's the relevant corresponding means
summary_percent$OutcomeMean <- c(0.0008868433,0.0008883278, 0.0003553007, 0.0004240056, 0.0003271193,0.0004071594,0.0003271193,0.0004071594,0.00006345647, 0.0002590487)
summary_percent$Mean <- (summary_percent$Mean / summary_percent$OutcomeMean) * 100
summary_percent$Lower <- (summary_percent$Lower / summary_percent$OutcomeMean) * 100
summary_percent$Upper <- (summary_percent$Upper / summary_percent$OutcomeMean) * 100
summary_percent$Test <- c("Placebo","Placebo","ENV","ENV","CY","CY","CSE","CSE","HR","HR")
#Row order
summary_percent$row_order <- interaction(summary_percent$Test, summary_percent$Area, sep = "_")
summary_percent$row_order <- factor(summary_percent$row_order,
                                    levels = c("Placebo_Agricultural", "Placebo_Forest",
                                               "ENV_Agricultural", "ENV_Forest",
                                               "CY_Agricultural", "CY_Forest",
                                               "CSE_Agricultural", "CSE_Forest",
                                               "HR_Agricultural", "HR_Forest")
)

# Custom labels to remove labels from x axis
x_labels <- c(
  " ", " ",
  " ", " ",
  " ", " ",
  " ", " ",
  " ", " "
)
#Rename the test names
section_labels <- data.frame(
  x = c("Placebo_Forest", "ENV_Forest", "CY_Forest", "CSE_Forest", "HR_Forest"),
  label = c("Placebo test", "Environmental", "Country year", "Country CSEs", "High risk")
)
#Define colors
effect_colors <- c("Forest" = "forestgreen","Agricultural" = "blue")
divider_pos <- c(2.5, 4.5, 6.5, 8.5)

##########
###Plot###
##########

setwd("~/OneDrive - Indiana University/FromGoogle/syd/plots/")
jpeg(filename = "sensitivityplot.jpeg", width = 8, height = 6, units = 'in', res = 500)
ggplot(summary_percent, aes(y = Mean, x = row_order, color = Area)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black", size = 1) +
  geom_vline(xintercept = divider_pos, linetype = "dotted", color = "grey40", size = 0.6) +
  scale_color_manual(values = effect_colors) +
  scale_x_discrete(labels = x_labels) +
  scale_y_continuous(
    limits = c(min(summary_percent$Lower) * 1.1, max(summary_percent$Upper) * 1.3),
    expand = c(0, 0),
    labels = scales::percent_format(scale = 1)
  ) +
  coord_flip(clip = "off") +
  # Section labels
  annotate("text",
           x =  c(1.5, 3.5, 5.5, 7.5, 9.5),
           y = -48,  #min(summary_percent$Lower) * 1.25,  # just outside plot on left
           label = section_labels$label,
           hjust = 0,
           size = 3.5,
           fontface = "bold.italic"
  ) +
  labs(y = "Effect size (%)", x = "", color = "Area type") +
  theme_minimal() +
  theme(
    panel.border = element_rect(colour = "black", fill = NA, size = 1),
    plot.title = element_text(hjust = 0.5)
  )

dev.off()





